body .main-container {
max-width: 100% !important;
width: 100% !important;
}
body {
max-width: 100% !important;
}
body, td {
font-size: 16px;
}
code.r{
font-size: 14px;
}
pre {
font-size: 14px
}Let us set some global options for all code chunks in this document.
knitr::opts_chunk$set(
message = FALSE, # Disable messages printed by R code chunks
warning = FALSE, # Disable warnings printed by R code chunks
echo = TRUE, # Show R code within code chunks in output
include = TRUE, # Include both R code and its results in output
eval = TRUE, # Evaluate R code chunks
cache = FALSE, # Enable caching of R code chunks for faster rendering
fig.align = "center",
out.width = "100%",
retina = 2,
error = TRUE
)
#klippy::klippy()
rm(list = ls())
set.seed(1982)Let us now load some required libraries.
# Load required libraries
library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)
library(plotly)
library(dplyr)
library(tidyr)
library(MASS)
library(glmnet)
library(car)
library(ggplot2)
library(plotly)
library(sf)
library(here)We now define some R functions with self-explanatory
names.
# Functions to manipulate the data
standardize <- function(vector) {return((vector - mean(vector)) / sd(vector))}
minmax_scaling <- function(x) {return((x - min(x)) / (max(x) - min(x)))}
convert_to_binary <- function(input_vector) {return(ifelse(input_vector != 0, 1, 0))}The following chunk uses the function
creates_covariates_on_mesh() defined on
Preprocessing/28Window.case.file4.addcovariatesonmesh.R to
build a mesh (on the same graph as the one used in this file) and
compute the value of all available covariates on the mesh nodes.
Go to Preprocessing 8.0.1 Raw
code to see the body of the
creates_covariates_on_mesh() function.
h = 0.05
source(here("Preprocessing/28Window.case.file4.addcovariatesonmesh.R"))
creates_covariates_on_mesh(h)
Sys.sleep(5)
################################################################################
################################# PREPARE THE DATA #############################
################################################################################
# loading the data
load(here("Data_files/data_on_graph_with_covariates_no_consecutive_zeros_partialtomtom.RData"))
load(here("Graph_objects/graph_construction_28_03_2024partialtomtomwhichlonglatsf.RData"))
load(here("Data_files/data_on_mesh_with_covariates_partialtomtom.RData"))
data = data_on_graph_with_covariates %>%
mutate(across(starts_with(c("class_", "upto")), list(ind = convert_to_binary))) %>%
mutate(across(c("bus", "signal", "stop", "crossing"), ~round(., 5))) %>%
mutate(across(c("SpeedLimit", "density_per_hour"), standardize))
mesh = data_on_mesh_with_covariates %>%
mutate(across(starts_with(c("class_", "upto")), list(ind = convert_to_binary))) %>% # this creates new columns
mutate(across(c("bus", "signal", "stop", "crossing"), ~round(., 5))) %>%
mutate(across(c("SpeedLimit", "density_per_hour"), standardize))Let us add the speed observations to the graph and build the mesh.
Observe that mesh already has the value of the available
covariates on the mesh nodes.
sf_graph$add_observations(data = data, group = "day", clear_obs = TRUE)
sf_graph$build_mesh(h = h)
Sys.sleep(5)stat.time.ini <- Sys.time()
################################################################################
################################# STATIONARY MODEL #############################
################################################################################
rspde_model_stat <- rspde.metric_graph(sf_graph,
parameterization = "matern",
nu = 0.5)
data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat,
repl = ".all",
loc_name = "loc")
cmp_stat = speed ~ -1 +
Intercept(1) +
SpeedLimit +
#density_per_hour +
#bus +
#signal +
#stop +
#crossing +
#upto1_ind +
#bus_number +
field(loc, model = rspde_model_stat,
replicate = data_rspde_bru_stat[["repl"]])
rspde_fit_stat <-
bru(cmp_stat,
data = data_rspde_bru_stat[["data"]],
family = "gaussian",
options = list(verbose = FALSE)
)
stat.time.fin <- Sys.time()
print(stat.time.fin - stat.time.ini)## Time difference of 35.79432 secsfit.rspde = rspde.result(rspde_fit_stat, "field", rspde_model_stat)
# summaries
# summary(rspde_fit_stat)
rspde_fit_stat$summary.fixedrbind(rspde_fit_stat$summary.hyperpar,
rspde_fit_stat$summary.hyperpar[1,]^(-0.5),
summary(fit.rspde))c(unlist(rspde_fit_stat$dic[c("dic", "dic.sat", "p.eff")]),
unlist(rspde_fit_stat$waic[c("waic", "p.eff")])) |> t() |> as.data.frame()rspde_fit_stat$mlik[2,1]## log marginal-likelihood (Gaussian)
## -55592.19nonstat.time.ini <- Sys.time()
################################################################################
############################# NON STATIONARY MODEL #############################
################################################################################
B.sigma = cbind(0, 1, 0, mesh$SpeedLimit, 0)
B.range = cbind(0, 0, 1, 0, mesh$SpeedLimit)
rspde_model_nonstat <- rspde.metric_graph(sf_graph,
start.theta = c(fit.rspde$summary.log.std.dev$mode,
fit.rspde$summary.log.range$mode,
rep(0, 2)),
B.sigma = B.sigma,
B.range = B.range,
parameterization = "matern",
nu = 0.5)
data_rspde_bru_nonstat <- graph_data_rspde(rspde_model_nonstat,
repl = ".all",
loc_name = "loc")
cmp_nonstat = speed ~ -1 +
Intercept(1) +
SpeedLimit +
#density_per_hour +
#bus +
#signal +
#stop +
#crossing +
#upto1_ind +
#bus_number +
field(loc, model = rspde_model_nonstat,
replicate = data_rspde_bru_nonstat[["repl"]])
rspde_fit_nonstat <-
bru(cmp_nonstat,
data = data_rspde_bru_nonstat[["data"]],
family = "gaussian",
options = list(verbose = FALSE)
)
nonstat.time.fin <- Sys.time()
print(nonstat.time.fin - nonstat.time.ini)## Time difference of 37.44091 secs# summaries
# summary(rspde_fit_nonstat)
rspde_fit_nonstat$summary.fixedrbind(rspde_fit_nonstat$summary.hyperpar,
rspde_fit_nonstat$summary.hyperpar[1,]^(-0.5),
summary(rspde.result(rspde_fit_nonstat, "field", rspde_model_nonstat)))c(unlist(rspde_fit_nonstat$dic[c("dic", "dic.sat", "p.eff")]),
unlist(rspde_fit_nonstat$waic[c("waic", "p.eff")])) |> t() |> as.data.frame()rspde_fit_nonstat$mlik[2,1]## log marginal-likelihood (Gaussian)
## -55598.8save.image(here(paste0("Models_output/", rmarkdown::metadata$title, ".RData")))